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BUCKLING  ANALYSIS  OF  CRACKED,  FLOATING  ICE  SHEETS 


M.D.  Adley  and  D.  S.  Sodhi 


INTRODUCTION 

Background 

The  design  load  for  offshore  structures  in  northern  waters  is  often  the 
horizontal  thrust  exerted  by  a  floating  ice  sheet.  The  maximum  ice  force 
that  may  be  developed  is  determined  by  either  the  environmental  driving 
forces  or  the  force  required  for  the  ice  sheet  to  fail.  The  smaller  of  these 
two  forces  will  govern.  Since  it  is  very  difficult  to  accurately  estimate 
the  magnitude  of  environmental  driving  forces  such  as  winds,  currents  and 
thermal  strains,  and  since  these  forces  may  generally  be  greater  than  the 
failure  load  of  the  ice  sheet,  the  force  required  for  the  ice  sheet  to  fail 
is  used  as  an  upper  limit  of  the  ice  forces. 

The  failure  load  of  the  ice  sheet  is  a  function  of  many  variables:  the 
strength  of  the  ice  sheet,  the  thickness  of  the  ice  sheet,  and  the  type  of 
structure-ice  interaction  are  some  of  the  important  factors.  These  factors 
will  determine  the  failure  mode  of  the  ice  sheet.  The  three  most  common 
modes  of  failure  for  the  ice  sheet  are  bending,  crushing  and  buckling.  The 
bending  mode  of  failure  is  generally  induced  by  a  fixed,  rigid  structure  with 
sloping  sides.  The  crushing  and  buckling  modes  of  failure  are  generally  in¬ 
duced  by  a  fixed,  rigid,  vertical  structure,  although  crushing  may  also  be 
caused  by  a  fixed,  flexible,  vertical  structure.  The  buckling  mode  of  fail¬ 
ure  is  most  likely  to  occur  when  the  aspect  ratio  (structure  width/ice  thick¬ 
ness)  is  large,  such  as  when  a  thin  ice  sheet  impinges  on  a  wide,  vertical 
structure.  Conversely  the  crushing  mode  of  failure  is  likely  to  occur  when 
the  aspect  ratio  is  small. 

When  a  floating  ice  sheet  is  loaded  with  an  in-plane  load,  cracks  may 
develop  from  the  point  of  loading  and  radiate  outwards  through  the  ice  sheet. 
As  a  result  of  those  cracks  the  domain  and  the  boundary  conditions  for  the 
differential  equation  governing  the  buckling  problem  are  altered.  Figure  1 


Structure 


a.  Symmetrical  configuration  of 
cracks  (  ot=  3)  . 


Structure 


b.  Unsymmetrical  configuration 
of  cracks  (  a*B) . 


Figure  1.  Plan  view  of  a  floating,  cracked  ice  sheet  and  a  structure  as  they 
push  each  other  to  produce  buckling  of  an  ice  sheet. 


shows  the  plan  view  of  a  floating  ice  sheet  that  is  pushed  against  a  struc¬ 
ture,  resulting  in  cracks  at  angles  a  and  3  from  the  edge  of  a  structure. 

The  orientation  of  cracks  will  be  divided  into  two  groups,  those  that  are 
symmetrical  about  the  line  of  loading  and  those  that  are  not  symmetrical 
about  the  line  of  loading.  Figures  la  and  b,  respectively,  show  examples  of 
a  symmetrical  and  an  unsymmetrical  configuration  of  cracks  relative  to  the 
structure. 

After  the  vertical  cracks  appear,  the  largest  force  the  ice  sheet  can 
exert  on  the  structure  is  the  failure  load  of  the  ice  sheet  still  in  contact 
with  the  structure.  If  the  aspect  ratio  is  large  enough  to  cause  buckling, 
the  largest  force  the  ice  sheet  can  exert  on  the  structure  is  the  buckling 
load  of  the  cracked,  floating  ice  sheets.  The  purpose  of  this  study  is  to 
determine  the  buckling  load  of  cracked,  floating  ice  sheets  as  opposed  to 
crushing  failure  of  ice  sheets.  The  results  of  theoretical  analyses  are  com¬ 
pared  to  those  determined  experimentally.  The  effect  of  the  ice  sheet  geome¬ 
try  on  the  buckling  load  is  determined  both  theoretically  and  experimentally. 
This  will  be  useful  in  determining  the  buckling  load  when  a  cracked  ice  sheet 
interacts  with  a  wide  structure. 

The  theoretical  part  of  the  analysis  consists  of  modeling  the  floating 
ice  sheet  as  a  thin,  homogeneous,  isotropic,  serai-infinite  plate  resting  on 
an  elastic  foundation.  The  force  exerted  by  the  elastic  foundation  is  as- 
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a=0°  0=0° 

Included  Angle  =  0° 


Included  Angle  =  30° 


Included  Angle=  120°  Included  Angle  =  150° 


Included  Angle  =90° 


Included  Angle  =  180° 


a.  Symmetrical  shapes  with  otf3  ^  90°.  b.  Symmetrical  shapes  with  otf 3>90°  . 


Figure  2.  Geometrical  shapes  of  the  ice  sheets  that  were  used  in  both  theo¬ 
retical  and  experimental  studies. 


sumed  to  be  linearly  proportional  to  the  deflection  of  the  ice  sheet,  and  the 
reactive  forces  at  adjacent  points  in  the  foundation  are  assumed  to  be  un¬ 
coupled,  i.e.  a  Winkler-type  foundation  is  used.  The  assumption  of  a  linear 
foundation  is  valid  as  long  as  the  ice  sheet  does  not  submerge  completely  un¬ 
der  the  water  surface  or  emerge  completely  out  of  the  water.  Thus,  this  as¬ 
sumption  is  valid  for  the  linear  stability  analysis  to  determine  the  bifurca¬ 
tion  load.  As  a  result  of  those  assumptions  a  differential  equation  can  be 
derived  that  describes  the  buckling  behavior  of  a  floating  ice  sheet.  This 
partial  differential  equation  is  solved  by  the  finite  element  method.  The 
solution  of  that  equation  yields  the  buckling  load  and  mode  of  buckling  of 
the  floating  ice  sheet. 

The  experimental  program  consisted  of  pushing  ice  sheets  of  different 
configurations  against  a  structure  and  measuring  the  buckling  load.  The  geo¬ 
metrical  shapes  considered  in  this  study  for  both  the  experimental  and  the 
theoretical  work  are  shown  in  Figure  2.  These  shapes  were  chosen  to  cover  a 
wide  range  of  angles;  the  results  for  the  intermediate  shapes  may  be  deter¬ 
mined  by  interpolation. 
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e.  Unsymmetrical  shapes  with  a=-30°  . 
Figure  2  (coat'd). 


Literature  review 

The  buckling  analysis  of  infinite  beams  on  elastic  foundations  was  pre¬ 
sented  by  Hetenyi  (1946).  Aside  from  this  early  work  the  buckling  analysis 
of  floating  ice  sheets  was  virtually  ignored  until  the  past  few  years. 

Sodhi  and  Hamza  (1977)  presented  a  stability  analysis  of  a  semi-infinite 
floating  ice  sheet,  which  was  performed  by  the  finite  element  method.  Kerr 
(1978)  presented  a  buckling  analysis  of  a  tapered  beam  of  ice.  Takagi  (1978) 
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solved  the  buckling  problem  of  an  infinite  floating  ice  sheet  uniformly 
stressed  along  the  periphery  of  an  internal  hole.  Wang  (1978)  presented  a 
buckling  anaylsis  of  a  semi-infinite  ice  sheet  moving  against  a  vertical, 
rigid,  cylindrical  structure;  the  finite  difference  method  was  used  for  that 
analysis. 

Nevel  (1979)  obtained  a  closed  form  solution  to  the  differential  equa¬ 
tion  that  governs  the  buckling  behavior  of  a  tapered  beam  of  ice  floating  on 
water.  This  analysis,  as  well  as  the  analysis  presented  by  Kerr  (1978),  was 
based  on  beam  theory.  Sodhi  (1979)  presented  a  buckling  analysis  of  a  semi¬ 
infinite  wedge-shaped  ice  sheet,  based  on  plate  theory. 

The  analysis  was  performed  by  the  finite  element  method.  In  that  report  Sod¬ 
hi  quantified  the  relationship  between  the  buckling  load  of  a  wedge-shaped 
ice  sheet  and  the  parameters  that  it  is  dependent  on,  i.e.  boundary  condi¬ 
tions,  aspect  ratio  and  included  angle. 

In  1980  a  review  of  buckling  analyses  was  presented  by  Sodhi  and  Nevel. 
This  report  compared  the  results  of  several  of  the  reports  previously  men¬ 
tioned. 

All  of  these  reports  were  based  on  elastic  linear  stability  theory. 
Therefore,  the  buckling  loads  computed  by  these  analyses  are  the  loads  at 
which  bifurcation  of  equilibrium  occurs,  not  the  failure  loads  of  the  ice 
sheets.  Kerr  (1981a)  presented  a  buckling  analysis  based  on  large  deforma¬ 
tion  theory  (i.e.  the  assumption  of  small  deformations  used  in  linear  stabil¬ 
ity  analyses  was  not  made)  . 

Sodhi  et  al.  (1982a)  conducted  an  experimental  study  of  the  buckling 
loads  required  for  ice  sheets  to  fail.  The  ice  sheets  were  pushed  against 
structures  of  different  widths,  and  the  force  required  for  the  ice  sheets  to 
fail  in  the  buckling  mode  was  measured.  Sodhi  (1983)  also  presented  a  dynam¬ 
ic  buckling  analysis  of  a  floating  ice  sheet. 

There  has  also  been  some  work  done  in  fields  other  than  ice  engineering 
that  addresses  the  problem  of  a  thin  plate  on  an  elastic  foundation.  Vlasov 
and  Leont’en  (1966)  discussed  buckling  analyses  of  beams,  plates  and  shells 
resting  on  elastic  foundations.  The  solutions  presented  in  that  paper  consi¬ 
dered  only  the  simply  supported  boundary  conditions.  The  results  of  a  few 
numerical  analyses  are  presented  in  the  Handbook  of  Structure  Stability  (Col¬ 
umn  Research  Committee  of  Japan  1971). 
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STATEMENT  OF  THE  PROBLEM 


The  governing  differential  equation  for  the  buckling  behavior  of  a  thin 
plate  on  an  elastic  foundation  can  be  written  as  follows  (Timoshenko  and  Gere 
1961): 

DV4w  +  Kw  =  Nxx|-f  +  2Nxy^-+Nyy|-f  (1) 

Bx  J  9x9y  9y 

where  w  *  out-of-plane  deflection  of  the  plate  (ice  sheet) 

x,  y  =  Cartesian  coordinates  on  the  mid-plane  of  the  plate  (ice 
sheet) 

V  =  biharmonic  operator 

^xx*  Nxy>  Nyy  *  in-plane  stress  resultants  (force  per  unit  length)  which  are 
linearly  dependent  on  the  total  in-plane  external  force  P. 

K  «  modulus  of  the  foundation  (specific  weight  of  water  in  the 
case  of  a  floating  ice  sheet) 

D  3  E  h3/l2(l-v2),  flexural  rigidity  of  the  plate  (ice  sheet) 

E  =  "effective"  elastic  modulus  of  the  ice  sheet 
h  =  ice  thickness 
v  ■  Poisson’s  ratio  for  ice. 

Equation  1  is  solved  by  first  solving  a  plane  stress  problem  to  deter¬ 
mine  the  in-plane  stress  resultants  ,  NXy  and  Nyy,  and  then  solving  the 
eigenvalue  problem  to  determine  the  buckling  loads  and  modes  of  buckling. 

The  nondimensional  form  of  the  buckling  load  is  found  by  normalizing  the  co¬ 
ordinates  x  and  y  with  respect  to  the  characteristic  length  of  the  plate  L 


L  =  L-- g-h-V-]lM 
12(1-vZ)K 


(2) 


to  obtain  jc  *  x/L  and  ^  *  y/L.  The  governing  equation  can  now  be  rewritten 
in  the  following  form: 


V  w  +  w 


~2 
<3  w 


BKL ' 


(Nxx  , 
-  dX 


+  2N 


»xy 


a  W 
<3x 


3 1 


(3) 


where  P  ■  total  in-plane  buckling  load 

B  =  width  over  which  P  is  distributed  at  the  boundary  of  the  domain 
Nxx  >NXy >Nyy  =  nondimensional  expressions  (Nj^B/P,  NxyB/P,  NyyB/P). 

Solutions  of  eq  3  that  satisfy  homogeneous  boundary  conditions  yield  the 
eigenvalue  P/BKL  (nondimensional  buckling  load)  and  the  corresponding  eigen¬ 
mode  (i.e.  mode  of  buckling) . 
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FINITE  ELEMENT  ANALYSIS 


The  stability  analysis  of  a  plate  on 
an  elastic  foundation  breaks  down  into 
two  separate  analyses,  an  in-plane  analy¬ 
sis  and  an  out-of-plane  analysis.  The 
in-plane  analysis  determines  the  stress 
distribution  within  the  plate,  which  is 
described  by  the  in- plane  stress  result¬ 
ants  obtained  in  the  analysis.  The  for¬ 
mulation  of  the  out-of-plane  problem  cul¬ 
minates  in  an  eigenvalue  problem.  The 
solution  of  this  problem  for  the  lowest 
eigenvalue  will  yield  the  lowest  buckling 
load  and  a  vector  of  displacement  repre 
senting  the  buckling  mode. 

In-plane  analysis 

The  element  used  in  the  in-plane  analysis  is  an  arbitrary  triangular 
element.  The  triangular  shape  was  chosen  for  its  effectiveness  in  modeling 
the  wedge-shaped  ice  sheets  that  will  be  analyzed.  The  element  and  its  coor¬ 
dinate  system  are  shown  in  Figure  3.  The  x,y  and  5, n  coordinate  systems  rep¬ 
resent  the  global  and  local  coordinate  systems,  respectively.  The  dimensions 
a,  b  and  c  shown  in  Figure  3  are  given  by 

a  =  [  (x  2~x  3 )  (x  2~x  i )  +  (y  2- y  3 )  (y  2-y  1 )  j/r 

b=  [  (x  3-x  i )  (x2~x  i )  +  (y3~y  iJ(y2~y  i)]/r  (4) 

c=  [(xrxi) Cy 3~y i )  -  (x3-xj(y2-yi)  ]/r 

where 


y 


Figure  3.  The  in-plane  and  out- 
of-plane  finite  element  and  its 
coordinate  system. 


r=  [ (x2~x i )2  +  (y2-yi)2]1^2  • 

The  x,y  coordinates  of  the  vertices  are  numbered  as  shown  in  Figure  3. 

The  strain  energy  expression  Ue  used  in  developing  the  plane  stress  ele¬ 
ment  may  be  written  as  follows  (Cowper  et  al.  1970): 


Ue  = 


Eh 


rr  rf-Su^  l(-3v'i2,  „  3u  3v  ,  1  N  ,-3u  ,  3v  ->2i .  /r.s 

J/  L(-s?J  +  1^)  +  2vT?^+  o  O-v)  Oar +-5?)  Jd5dn  (5) 


2(i-v2)  JJ  LV95;  ’  3C  3n  '  2  Vi  v/  ^3n  ’  35 


where  u  is  the  displacement  in  the  5  direction,  and  v  is  the  displacement  in 
the  n  direction. 
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The  displacement  functions  used  to  represent  the  u  and  v  displacements 
within  an  element  are 


The  finite  element  derived  using  these  displacement  functions  is  a  20- 
degree-of-freedom  element.  There  are  six  degrees  of  freedom  at  each  of  the 


are  also  two  degrees  of  freedom  (u  and  v)  at  the  centroid  of  the  element, 
which  are  eliminated  by  a  static  condensation  procedure. 

The  displacement  functions  used  in  developing  this  element  will  lead  to 


strain  energy  convergence  rates  proportional  to  n~6,  where  n  is  the  number  of 


elements  per  side  of  a  structure. 

By  utilizing  these  shape  functions  and  the  strain  energy  expression  pre¬ 
sented  above,  the  elements  of  the  stiffness  matrix  can  be  represented  by 
closed  form  expressions.  For  a  complete  derivation  of  the  element,  see  Cow- 
per  et  al.  (1970) . 

Out-of-plane  analysis 

The  element  used  in  the  out-of-plane  analysis  is  also  an  arbitrary  tri¬ 
angular  element.  The  coordinate  system  for  the  out-of-plane  element  is  iden¬ 
tical  to  that  used  for  the  in-plane  element  (Fig.  3). 

The  strain  energy  expression  for  the  stability  analysis  may  be  written 
as  follows  (Gallagher  1975): 


(7) 


where  w  is  the  out-of-plane  deflection  and  a  refers  to  the  stresses  resulting 
from  the  applied  axial  load  (tensile  stress  is  defined  as  positive). 

The  displacement  function  used  to  represent  the  out-of-plane  deflection 
within  an  element  is 


w(5,t))  =  ai  +  a2^  +  a  3TI  +  at+?2  +  asSn  +  aen2  +  a7?3  +  a8^2n  +  a  g?n2 

+  aion3  +  an?4  +  ai253Ti  +  ai3?2n2  +  am?n3  +  aisn4  +  ai6?5  (8) 

+  ai7?3ri2  +  ai8?2Ti3  +  ai9?rj4  +  a20r'5* 
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By  substituting  eq  8  into  eq  7,  we  may  derive  the  stiffness  and  geometric 
matrices  for  the  element. 

The  finite  element  derived  using  the  displacement  function  presented  in 
eq  8  is  an  18-degree-of-freedom  element.  There  are  six  degrees  of  freedom  at 
each  vertex  of  the  triangle.  These  are  the  transverse  deflection  and  its 
first  and  second  derivatives,  i.e. 

3w  3w  32w  32w  32w 

w  H  3n  "a?  3£3n  * 

The  displacement  function  used  in  developing  this  element  is  a  quintic  poly¬ 
nomial  in  x  and  y.  Therefore,  the  deflection  along  any  edge  of  the  element 
varies  as  a  quintic  polynomial  in  the  edgewise  coordinate. 

The  general  quintic  polynomial  in  two  variables  depends  on  21  constants. 
Since  there  are  only  18  degrees  of  freedom,  3  additional  constraints  may  be 
satisfied:  the  slope  normal  to  each  edge  must  be  a  cubic  function  of  the 
edgewise  coordinate.  The  four  coefficients  of  the  cubic  polynomial  are 
uniquely  determined  by  the  slope  normal  to  the  edge  and,  by  the  twist,  at 
each  of  the  two  terminal  vertices.  Therefore,  the  continuity  of  displace¬ 
ments  and  normal  slopes  are  maintained  across  interelement  boundaries. 

The  accuracy  of  the  out-of-plane  element  is  comparable  to  the  accuracy 
of  the  in-plane  element.  This  is  very  desirable  because  the  accuracy  of  the 
results  of  the  in-plane  analysis  will  affect  the  accuracy  of  the  out-of-plane 
analysis.  For  a  complete  derivation  of  the  element,  see  Cowper  et  al. 

(1969). 

The  finite  element  formulation  of  the  linear  stability  analysis  of  a 
plate  on  an  elastic  foundation  may  be  derived  by  utilizing  the  theory  of  sta¬ 
tionary  potential  energy.  By  employing  the  principle  of  stationary  potential 
energy  and  following  a  well-established  procedure  (Gallagher  1973),  the  fol¬ 
lowing  equation  may  be  derived: 

[K]{Af}=  A[KG]{Af}  (9) 

where 

[K]  *  stiffness  matrix 

[KG]  =  geometric  stiffness  matrix 

{ Af }  =  vector  of  displacements  (the  eigenvector) 

X  =  multiplication  factor  (the  eigenvalue)  . 
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Equation  9,  which  is  formulated  in  the  out-of-plane  analysis,  is  in  the 
form  of  an  eigenvalue  problem.  The  solution  of  eq  9  yields  the  buckling  load 
X  and  the  buckling  mode  {Af}. 

Two  finite  element  computer  programs  have  been  written  for  this  project. 
The  first  program  (P0EF1)  was  written  to  solve  the  in-plane  stress  problem. 
P0EF1  utilizes  the  in-plane  stress  element  previously  discussed.  The  second 
program  (P0EF2)  was  written  to  formulate  the  out-of-plane  problem.  P0EF2 
utilizes  the  out-of-plane  element  previously  discussed.  A  third  program  was 
an  in-plane  stress  program  that  utilized  a  triangular  element  with  a  linear 
displacement  function.  This  program  was  taken  from  Desai  and  Abel  (1972). 

The  finite  elements  used  in  P0EF1  and  P0EF2  were  chosen  for  their  proven 
high  accuracy  and  good  error-convergence  rates.  Since  the  out-of-plane  anal¬ 
ysis  is  affected  by  any  errors  in  the  in-plane  analysis,  it  is  also  desirable 
that  the  elements  used  both  have  the  same  strain-energy  convergence  rates. 

The  finite  element  programs  used  for  this  project  were  thoroughly  test¬ 
ed.  Each  of  the  three  programs  was  used  to  analyze  problems  for  which  an  an¬ 
alytical  solution  was  available.  The  results  were  very  good.  The  solutions 
given  by  the  computer  programs  agreed  with  the  analytical  solutions,  and  the 
convergence  properties  appeared  to  be  very  good. 

Reduction  of  the  eigenvalue  problem 

The  out-of-plane  analysis  culminates  in  an  eigenvalue  problem  in  the 
form  of  eq  9 ,  which  is  solved  by  two  methods.  The  first  method  involves 
solving  the  full-sized  eigenvalue  problem.  This  requires  a  very  effective 
solution  algorithm  and  also  places  restrictions  on  the  size  of  the  problems. 
Therefore,  only  selected  problems  are  solved  by  the  first  method. 

The  second  method  includes  an  algorithm  to  reduce  the  size  of  the  eigen¬ 
value  problem.  In  this  method,  known  as  the  Guyan  reduction  method  (Guyan 
1965),  some  degrees  of  freedom  are  designated  as  masters  and  some  as  slaves. 
The  degrees  of  freedom  retained  in  the  reduced  eigenvalue  problem  are  defined 
as  the  masters.  The  Guyan  reduction  method  utilizes  the  static  relationship 
between  the  masters  and  the  slaves  to  remove  the  slave  degrees  of  freedom 
while  preserving  the  total  strain  energy  of  the  structure.  The  problem  of 
selecting  the  masters  is  overcome  by  employing  a  recently  developed  algorithm 
that  allows  an  analytical  selection  of  the  masters.  The  second  method  signi¬ 
ficantly  reduces  the  size  of  the  eigenvalue  problem  while  preserving  the  ac¬ 
curacy  of  the  eigenvalues  of  interest. 
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The  following  is  a  presentation  of  the  algorithm  used  for  the  analytical 
selection  of  the  masters.  This  algorithm,  along  with  a  more  thorough  expla¬ 
nation  of  its  use,  can  be  found  in  Shah  and  Raymund  (1982).  The  eigenvalue 
problem  to  be  solved  is  in  the  form  of  eq  9. 

Step  1.  Select  a  cut-off  value  of  A,  to  be  known  as  Ac. 

Step  2.  Find  the  degree  of  freedom  for  which  the  ratio  of  ^±±/^G±±  =  X 

is  the  largest.  If  several  degrees  of  freedom  have  the  same  value  of  the 
ratio,  choose  the  one  with  the  lowest  index. 

Step  3.  If  this  value  is  greater  than  Ac,  condense  this  degree  of  free¬ 
dom  from  the  stiffness  and  geometric  stiffness  matrices  by  the  Guyan  reduc¬ 
tion  procedure. 

Step  4.  Apply  steps  2  and  3  to  the  reduced  matrices  obtained. 

Step  5.  If  the  greatest  value  of  the  ratio  K^-j/KG-q  is  _<  Ac,  stop  the 

procedure  at  this  point  and  use  the  reduced  matrices  obtained  to  calculate 
the  eigenvalues  and  eigenvectors  of  interest. 

This  algorithm  has  been  used  in  dynamic  natural  frequency  analyses  and 
has  met  with  a  great  deal  of  success  (Shah  and  Raymund  1982).  P0EF2  utilizes 
this  algorithm  to  reduce  the  size  of  the  eigenvalue  problem. 

Development  of  the  finite  element  model 

The  selection  of  the  finite  element  model  is  a  crucial  step  in  an  analy¬ 
sis  because  the  results  can  be  greatly  affected  by  the  characteristics  of  the 
model.  The  finite  element  models  used  for  this  project  were  developed  by  us¬ 
ing  several  different  models  to  analyze  the  same  problem.  The  error  intro¬ 
duced  by  varying  the  level  of  discretization  was  determined  by  comparing  the 
results  from  the  different  models.  The  models  used  in  the  analysis  were 
chosen  for  their  ability  to  accurately  describe  the  behavior  of  the  ice  sheet 
in  spite  of  the  following  considerations. 

1)  A  finite-sized  model  is  used  to  represent  the  behavior  of  a  semi¬ 
infinite  ice  sheet.  Therefore,  the  model  must  be  large  enough  so  that  the 
solution  is  not  significantly  affected  by  changes  in  boundary  conditions  at 
the  edges,  which  are  supposed  to  be  at  infinity.  However,  the  model  has  to 
be  small  enough  so  that  the  solution  will  be  computationally  feasible. 

2)  The  discretization  of  the  structure  must  be  designed  so  that  areas 
with  large  stress  and  displacement  gradients  are  highly  discretized. 

3)  The  in-plane  stress  programs  calculate  the  in-plane  stress  resultants 
at  the  centroids  of  the  finite  elements.  Therefore,  if  the  element  is  too 


11 


large,  the  values  of  the  in-plane  stress  resultants  will  not  be  representa¬ 
tive  of  the  stresses  in  the  area  covered  by  the  element. 

It  is  obvious  that  the  solution  is  affected  by  the  size  and  number  of 
elements  in  the  model.  It  is  advantageous,  then,  to  analyze  the  symmetrical¬ 
ly  shaped  ice  sheets  and  the  unsymmetrically  shaped  ice  sheets  separately. 

Finite  element  model  for  the  symmetrical  shapes 

Symmetrically  shaped  ice  sheets,  i.e.  wedges,  have  been  analyzed  previ¬ 
ously  by  other  researchers.  However,  the  previous  analyses  used  closed  form 
expressions  to  describe  the  distribution  of  in-plane  stresses.  In  the  analy¬ 
sis  presented  here,  a  separate  in— plane  finite  element  analysis  is  performed 
to  ascertain  the  distribution  of  in-plane  stresses  within  the  element.  Be¬ 
cause  of  the  generality  of  the  finite  element  method,  loading  distributions 
can  be  chosen  that  more  closely  simulate  conditions  that  might  occur  in  na¬ 
ture  or,  as  in  this  case,  conditions  existing  during  model  studies. 

The  finite  element  model  used  to  analyze  one  of  the  symmetrical  shapes 
is  shown  in  Figure  4.  Figure  5  shows  the  boundary  conditions  and  the  load 
distribution  used  in  the  in-plane  analysis.  The  edge  that  is  supposed  to  go 
to  infinity  is  fixed  in  the  x  direction  but  is  free  to  displace  in  the  y  dir¬ 
ection.  Since  the  wedge  is  symmetrical  about  its  longitudinal  axis,  only 
half  of  the  wedge  is  modeled.  Along  the  centerline  the  model  is  fixed  in  the 
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Figure  5.  Boundary  conditions  and 
load  distribution  used  in  the  in¬ 
plane  analysis  of  the  symmetrical 
shapes.  The  rollers  represent  zero 
normal  displacement  and  zero  tangen¬ 
tial  force  at  the  boundary.  The 
uniformly  distributed  in-plane  load 
intensity  P/B  is  applied  at  the  ice/ 
structure  interface. 


Figure  6.  Boundary  conditions  used 
in  the  out-of-plane  analysis  of  the 
symmetrical  shapes.  The  free  bound¬ 
ary  condition  signifies  zero  bending 
moment  and  zero  effective  shear  force 
at  the  boundary.  The  fixed  boundary 
condition  signifies  zero  deflection 
and  zero  normal  slope. 


y  direction  due  to  symmetry  but  is  free  to  deform  along  the  x  axis.  The 

boundary  conditions  on  the  loaded  edge  and  the  tapered  edge  are  modeled  as 
free. 

The  loading  distribution  used  in  the  analysis  is  uniform.  Since  the  re¬ 
sults  of  the  theoretical  analysis  will  be  compared  to  the  experimental  re¬ 
sults  from  the  model  study,  the  loading  was  chosen  to  represent  the  loading 
that  the  ice  sheets  experienced  during  the  tests.  The  true  loading  distribu¬ 
tion  the  ice  sheets  experienced  is  actually  unknown,  but  the  authors  feel 
that  the  uniformly  distributed  load  is  a  good  assumption. 

There  are  three  boundary  conditions  used  in  the  out-of-plane  analysis 
(Fig.  6): 

Free:  moment  (M)  =  0,  shear(V)  =  0 
Fixed:  displacement  (w)  =  0,  slope  (8w/9x)  =0 

Hinged:  displacement  (w)  =  0,  moment  (M)  =  0  . 

The  edge  that  is  supposed  to  go  to  infinity  is  modeled  as  fixed.  Since  the 
ice  sheet  is  symmetrical  about  its  longitudinal  axis,  only  half  of  the  ice 
sheet  is  analyzed.  The  boundary  condition  along  the  centerline  is  free. 
However,  due  to  symmetry  the  slope  normal  to  the  centerline  is  set  to  zero. 
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The  boundary  condition  along  the  tapered  edge  is  modeled  as  free.  The  bound¬ 
ary  along  the  loaded  edge  is  modeled  as  a  hinged  boundary. 

Finite  element  model  for  the  unsymmetrical  shapes 

A  finite  element  model  used  to  analyze  one  of  the  unsymmetrical  shapes 
is  shown  in  Figure  7.  The  boundary  conditions  and  load  distribution  used  in 
the  in-plane  analysis  of  an  unsymmetrical  ice  sheet  are  shown  in  Figure  8. 

The  edge  that  is  supposed  to  go  to  infinity  is  modeled  as  fixed  in  the  x 
direction,  but  it  is  free  to  displace  in  the  y  direction.  The  loaded  edge  is 
modeled  as  free  to  move  in  both  the  x  and  y  directions.  The  in-plane  dis¬ 
placements  in  the  direction  perpendicular  to  the  cracks  were  assumed  to  be 
zero,  as  the  ice  across  the  crack  can  prevent  the  displacement  in  the  normal 
but  not  in  the  tangential  direction. 

The  shape  of  the  unsymmetrical  ice  sheet  and  the  loading  configuration 
are  such  that  the  ice  sheet  will  experience  a  moment  as  well  as  a  compressive 
load  (Fig.  9).  The  in-plane  boundary  condition  at  the  cracks  is  to  allow 
displacement  in  the  tangential  direction  (slip)  and  no  displacement  in  the 
normal  direction  due  to  the  presence  of  adjoining  ice.  The  magnitude  of  the 
moment  experienced  by  the  ice  sheet  depends  on  the  shape  of  the  distributed 
loading  applied.  The  authors  chose  to  model  this  loading  as  a  trapezoidal 
load.  There  are  now  two  unknowns  that  pertain  to  the  loading  distribution: 
the  magnitude  of  the  uniform  portion  of  the  trapezoidal  load  and  the  slope  of 
the  top  side  of  the  trapezoid  (Fig.  8).  These  two  variables  (b  and  0)  cannot 
be  solved  for  without  additional  information. 


Figure  7.  One  of  the  finite  element  models  used  to  analyze 
the  unsymmetrical  shapes. 
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Figure  8,  Boundary  conditions  and  in- plane 
loading  distributions  used  in  the  in-plane 
analysis  of  the  unsymmetr ical  shapes. 


Figure  9.  Moment  induced  by  the  unsymmetrical 
shape  and  the  partially  confined  in-plane  bound¬ 
ary  conditions  used  as  a  consequence  of  that 
moment. 


Two  methods  of  dealing  with  this  problem  were  investigated.  In  one  case 
the  additional  information  needed  to  solve  for  b  and  0  was  taken  from  the  ex¬ 
perimental  force  records.  In  the  second  method  the  loading  distributions 
that  gave  the  largest  and  smallest  buckling  pressures  were  used.  This  method 
yielded  two  curves,  one  representing  the  largest  force  that  the  ice  sheet 
would  exert  on  the  structure  and  the  other  the  minimum  force  required  for  the 
ice  sheet  to  fail.  Results  were  computed  for  both  of  these  methods.  How- 
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Hinged  Figure  10.  Boundary  conditions  used 

for  the  out-of-plane  analysis  of  the 
— unsymmetrical  shapes. 

ever,  the  authors  favor  the  second  method  since  more  useful  information  can 
be  presented  without  introducing  further  sources  of  error. 

An  example  of  the  boundary  conditions  used  in  the  out-of-plane  analysis 
for  the  unsymmetrical  shapes  is  shown  in  Figure  10.  The  boundary  conditions 
used  for  the  out-of-plane  analysis  are  the  following:  hinged  at  the  loaded 
edge,  fixed  at  the  edge  that  is  supposed  to  go  to  infinity,  and  free  on  the 
remaining  two  edges. 

EXPERIMENTAL  PROCEDURE 

The  model  study  presented  in  this  report  was  performed  in  the  refriger¬ 
ated  test  basin  at  CRREL  (Sodhi  and  Adley  1984).  The  ice  sheets  tested  were 
grown  by  seeding  and  freezing  a  solution  of  1%  urea  in  water  at  an  ambient 
temperature  of  -12°C.  Prior  to  the  tests  the  refrigerated  room  was  warmed  up 
and  the  temperature  of  the  ice  was  allowed  to  stabilize  at  0°C.  The  ice 
tested  can  be  characterized  as  columnar  ice  with  horizontal  c-axes. 

The  model  study  consisted  of  monitoring  the  horizontal  forces  while 
either  pushing  a  structure  against  a  stationary  ice  sheet  or  pushing  an  ice 
sheet  against  a  stationary  structure.  The  relative  velocity  between  the 
structure  and  the  ice  sheet  was  held  constant  at  1  cm/s.  In  both  cases  the 
length  of  the  structur e/ice  interface  was  1.85  m.  The  experimental  set-up  is 
shown  in  Figure  11. 

Prior  to  the  tests  the  thickness,  flexural  strength  and  characteristic 
length  of  the  ice  sheet  were  measured.  The  characteristic  length  was  deter¬ 
mined  so  that  the  experimental  results  could  be  compared  to  the  theoretical 
results.  The  procedure  involves  placing  dead  weights  in  discrete  increments 
on  the  ice  sheet  and  monitoring  its  deflection.  This  procedure  is  described 
by  Sodhi  et  al.  (1982a). 
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a.  Side  view. 


Carriage 


let  Shtel 


Instrumented 
Support  Bar 


b.  Schematic  view. 


Figure  11.  Experimental  set-up  used  in  the  model  study. 


The  edge  of  the  ice  sheet  was  sawed  parallel  to  the  structure  prior  to 
testing  to  ensure  uniform,  continuous  contact.  The  geometric  shape  of  the 
ice  sheet  was  modified  by  cutting  two  slots  with  a  saw.  An  example  of  these 
slots  is  shown  in  Figure  12  (the  slots  are  marked  by  the 
black-and-white-striped  poles) •  The  far  edge  of  the  ice  sheet  was  left 
intact.  A  3-mm-thick  rubber  pad  was  attached  to  the  structure  where  the 
structure/ice  interface  would  occur  during  the  tests.  The  pad  prevented  the 
edge  of  the  ice  sheet  from  moving  up  or  down  but  left  it  free  to  rotate, 
creating  a  hinged  boundary  condition.  This  condition  was  chosen  because  it 
could  be  simulated  in  the  experiments  quite  accurately. 

The  only  difference  between  tests  was  the  orientation  of  the  two  slots, 
which  were  cut  to  modify  the  shape  of  the  ice  sheet.  Thus,  the  only 
parameter  varied  was  the  geometric  shape.  The  buckling  load  was  assumed  to 
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a.  Symmetrical  ice  sheet  prior  to  testing. 


b.  Symmetrical  ice  sheet  during  testing. 


Figure  12.  Ice  sheet  buckling  tests. 
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c.  Un symmetrical  ice  sheet  during  testing. 


d.  Unsymmetrical  ice  sheet  after  testing. 
Figure  12  (contfd). 
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Figure  13*  Typical  ice  forces  measured  at 
the  two  supports  in  the  model  study. 


be  the  largest  (peak)  force  exerted  on  the  structure  by  the  ice  sheet  as  de¬ 
termined  from  the  experimental  force  records  (Fig.  13). 

RESULTS 

2 

The  nondimensional  buckling  load  (P/BKL  )  is  a  function  of  the  boundary 
conditions,  the  included  angle,  and  the  aspect  ratio.  For  the  theoretical 
analysis  the  boundary  conditions  and  the  aspect  ratio  were  held  constant. 
Thus,  the  results  presented  will  show  the  effect  that  the  geometry  of  the  ice 
sheet  has  on  the  nondimensional  buckling  load.  The  boundary  conditions  used 
in  the  theoretical  analysis  were  chosen  to  simulate  the  boundary  conditions 
used  in  the  model  study,  so  the  results  of  the  theoretical  analysis  can  be 
compared  to  the  experimental  results. 

Symmetrical  shapes 

The  results  of  the  analyses  of  the  symmetrically  shaped  ice  sheets  are 
presented  in  Figure  14.  The  experimental  results  are  also  presented  in  Table 
1,  along  with  data  on  ice  sheet  and  test  conditions.  In  approximately  90%  of 
the  tests  the  aspect  ratio  of  the  ice  sheet  is  very  close  to  that  used  in  the 
theoretical  analysis  (4.3).  Thus,  there  is  a  strong  basis  for  comparison  be¬ 
tween  most  of  the  points  and  the  curves  developed  in  the  theoretical  analy¬ 
sis. 

The  three  curves  presented  in  Figure  14  represent  separate  theoretical 
analyses  but  use  the  same  aspect  ratio  and  boundary  conditions.  Two  of  the 
curves,  representing  the  results  of  two  finite  element  stability  analyses, 
were  computed  for  this  project.  For  the  first  analysis  (represented  by  curve 
1)  a  triangular-shaped  element  utilizing  a  cubic  displacement  function  was 
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Table  1 


Test 

no. 


14 

15 

16 

17 

18 

19 

20 
21 

25 

26 

27 

28 

29 

30 

31 

32 


Results  of  the  experiment  on  symmetrically  shaped  ice  sheets. 


Included 


Angle 

(degrees) 

h 

(cm) 

L 

(m) 

B 

(m) 

180 

3.4 

0.46 

1.85 

0 

3.4 

0.46 

1.85 

60 

3.4 

0.46 

1.85 

90 

3.4 

0.39 

1.85 

120 

3.4 

0.39 

1.85 

30 

3.4 

0.39 

1.85 

0 

3.4 

0.35 

1.85 

180 

3.4 

0.35 

1.85 

0 

2.8 

0.475 

1.85 

180 

2.8 

0.475 

1.85 

60 

2.8 

0.45 

1.85 

90 

2.8 

0.435 

1.85 

120 

2.8 

0.43 

1.85 

150 

2.8 

0.43 

1.85 

30 

2.8 

0.43 

1.85 

0 

2.8 

0.43 

1.85 

P 

(N) 

B/L 

P 

- 2 

BKL 

10,353.0 

4.02 

2.697 

9,468.0 

4.02 

2.466 

8,891.0 

4.02 

2.316 

8,476.0 

4.74 

3.072 

9,267.0 

4.74 

3.358 
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Figure  14.  Results  of  the  theoretical 
(lines)  and  experimental  (points)  work 
on  symmetrical  shapes.  The  assumptions 
for  these  analyses  are  described  in  the 
text. 
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Figure  15.  Plot  of  nondime ns ional  buckling 
pressure  (P/BKL  )  vs  included  angle  at- 3  for 
different  boundary  conditions  at  the  ice/ 
structure  interface  and  for  different  aspect 
ratios  (B/L). 

employed  for  the  in-plane  stress  analysis*  This  analysis  used  the  computer 
programs  written  for  this  project  (P0EF1  and  P0EF2).  In  the  second  analysis 
(represented  by  curve  2)  a  triangular-shaped  element  utilizing  a  linear  dis¬ 
placement  function  was  used  for  the  in-plane  stress  analysis*  This  analysis 
used  a  plane  stress  finite  element  program  taken  from  Desai  and  Abel  (1972) 
and  the  out-of-plane  program  written  for  this  project  (P0EF2) •  Both  of  these 
analyses  used  a  computer  program  presented  by  Gladwell  and  Tahbildar  (1972) 
that  solves  the  algebraic  eigenvalue  problem* 

The  third  curve  (represented  by  curve  3)  was  computed  by  Sodhi  (1979). 
Sodhi  used  a  Fourier  decomposition  and  the  finite  element  method  to  solve  the 
buckling  problem  of  the  symmetrical  shapes.  Sodhi’ s  analysis  assumes  a  radi¬ 
al  stress  field  for  the  in-plane  stress  distribution. 

For  a  symmetrical  configuration  of  cracks  the  results  of  the  buckling 
analysis  of  the  wedge-shaped  ice  sheet  (Sodhi  1979)  have  been  plotted  in  Fig- 
ure  15,  which  shows  the  normalized  buckling  loads  (P/BKL  )  with  respect  to 
the  wedge  angle  (otf  3)  for  different  aspect  ratios  and  boundary  conditions  at 
the  ice/structure  interface. 

Unsymmetrical  shapes 

The  results  of  the  analyses  of  the  unsymmetrical ly  shaped  ice  sheets  are 
presented  in  Figure  16.  The  test  results  are  presented  in  Table  2. 
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Table  2.  Results  of  the  experiment  on  unsymmetrically  shaped  ice  sheets. 
Included 


Test 

no. 

angles 

(degrees) 

a  8 

h 

(cm) 

L 

(m) 

B 

(o) 

P 

(N) 

B 

L 

P 

BKL* 

33 

0 

0 

3.03 

0.35 

1.85 

4,279.0 

5.29 

1.925 

34 

0 

90 

3.03 

0.35 

1.85 

5,034.0 

5.29 

2.265 

35 

0 

30 

3.03 

0.35 

1.85 

4,168.0 

5.29 

1.875 

36 

0 

60 

3.03 

0.35 

1.85 

4,579.0 

5.29 

2.060 

37 

0 

45 

3.03 

0.30 

1.85 

4,742.0 

6.17 

2.904 

38 

-30 

60 

3.03 

0.30 

1.85 

4,775.0 

6.17 

2.925 

55 

-15 

15 

3.10 

0.39 

1.85 

8,592.0 

4.74 

3.114 

56 

-15 

30 

3.10 

0.39 

1.85 

7,911.0 

4.74 

2.867 

57 

-15 

60 

3.10 

0.39 

1.85 

6,305.0 

4.74 

2.285 

58 

-15 

90 

3.10 

0.37 

1.85 

6,262.0 

5.00 

2.521 

59 

-30 

30 

3.10 

0.35 

1.85 

3,614.0 

5.29 

1.626 

60 

-30 

60 

3.10 

0.35 

1.85 

3,658.0 

5.29 

1.646 

61 

-30 

90 

3.10 

0.35 

1.85 

3,405.0 

5.29 

1.532 

69 

0 

0 

2.87 

0.32 

1.85 

5,097.0 

5.78 

2.744 

71 

0 

60 

2.87 

0.32 

1.85 

6,671.0 

5.78 

3.591 

72 

0 

90 

2.87 

0.32 

1.85 

4,245.0 

5.78 

2.285 

The  curves  presented  in  Figures  16a  and  b  are  the  results  of  the  theo¬ 
retical  analyses.  The  same  aspect  ratio  and  boundary  conditions  were  used 
for  each  analysis,  and  the  aspect  ratio  was  the  same  as  that  used  for  the 
symmetrical  analyses. 

The  curves  were  computed  using  an  in-plane  finite  element  program  pre¬ 
sented  by  Desai  and  Abel  (1972)  and  the  out-of-plane  program  written  for  this 
project  (P0EF2).  The  program  presented  by  Desai  and  Abel  (1972)  utilizes  a 
triangular-shaped  element  derived  by  using  a  linear  displacement  function. 

The  eigenvalue  problem  was  solved  by  using  the  computer  program  presented  by 
Gladwell  and  Tahbildar  (1972). 

The  three  curves  presented  in  Figure  16  were  computed  using  different 
loading  distributions.  Curves  1  and  2  represent,  respectively,  the  largest 
load  the  ice  sheet  can  withstand  prior  to  buckling  and  the  minimum  load  re¬ 
quired  to  buckle  the  ice  sheet.  Curve  3  is  computed  by  applying  the  loading 
distributions  obtained  from  experimental  force  records. 

There  are  no  results  available  for  the  unsymmetrical  shapes  when  a  = 

-30°  due  to  numerical  difficulties  in  obtaining  a  solution.  The  cause  for 
these  difficulties  is  not  known. 
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0 


a.  a  =  0°  . 


b.  a  =  -15. 


c. 


a  -  -30° . 


Figure  16.  Results  of  the  theoretical  (lines)  and 
experimental  (points)  work  on  unsymmetrical  shapes. 
The  theoretical  results  were  obtained  by  finite  ele¬ 
ment  analysis  using  different  assumptions  for  the 
distribution  of  in-plane  load. 
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Mode  shapes 

Since  the  load  was  applied  slowly  to  the  ice  sheets  tested  in  the  model 
study  (1  cm/s),  all  of  the  ice  sheets  buckled  in  the  first  mode.  The  shape 
of  the  first  mode  of  buckling  predicted  by  the  theoretical  analysis  appears 
to  agree  well  with  the  mode  shapes  observed  in  the  experimental  work. 

DISCUSSION  AND  CONCLUSIONS 

The  theoretical  analyses  presented  in  this  report  were  done  for  one  val¬ 
ue  of  the  aspect  ratio  (structure  width/characteristic  length).  However,  the 
close  agreement  between  the  results  computed  for  this  project  and  the  results 
computed  by  Sodhi  (1979)  and  between  the  theoretical  and  experimental  results 
of  this  project  indicate  that  the  analysis  developed  could  be  used  success¬ 
fully  for  ice  sheets  with  different  aspect  ratios.  The  comparison  between 
the  results  presented  in  this  paper  and  Sodhi' s  (1979)  results  also  indicates 
that  the  assumption  of  the  radial  stress  field  for  the  distribution  of  the 
in-plane  stresses  is  valid  for  the  symmetrical  shapes. 

The  authors  believe,  after  comparing  the  results  of  the  analyses  of  the 
symmetrical  shapes,  that  the  use  of  the  higher-order  plane  stress  finite  ele¬ 
ment  was  not  warranted.  The  results  of  the  stability  analyses  were  not  ad¬ 
versely  affected  by  using  the  less-accurate  plane  stress  computer  program 
taken  from  Desai  and  Abel  (1972),  so  it  was  used  for  the  analyses  of  the  un- 
symmetrical  shapes  because  it  is  less  expensive  to  run. 

The  results  of  the  theoretical  and  experimental  analyses  show  that  the 
nondimensional  buckling  load  for  a  symmetrical  configuration  increases  as  the 
included  angle  increases.  Thus,  the  configuration  with  the  largest  buckling 
load  is  the  symmetrical  shape  with  an  included  angle  of  180° . 

The  in-plane  analysis  of  the  unsymmetrical  shapes  became  quite  complex 
because  of  1)  the  moment  induced  by  the  unsymmetrical  geometry,  2)  the  abili¬ 
ty  of  the  cracks  to  transmit  compressive  normal  forces  created  by  that  mo¬ 
ment,  and  3)  the  fact  that  the  loading  distribution,  and  therefore  the  magni¬ 
tude  of  the  applied  normal  force  and  moment,  is  an  unknown.  Because  of  that 
moment,  the  in— plane  analysis  was  done  by  modeling  the  unsymmetrical  shapes 
as  partially  confined;  i.e.  no  displacements  were  allowed  perpendicular  to 
one  longitudinal  edge. 

Since  the  actual  loading  distribution  experienced  by  the  ice  sheet  is  an 
unknown,  the  analysis  used  two  assumed  loading  distributions  and  one  loading 
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distribution  determined  from  the  experimental  force  records.  The  two  assumed 
loading  distributions  were  chosen  to  give  the  maximum  and  minimum  buckling 
loads.  The  loading  distribution  determined  from  the  force  record  approxi¬ 
mates  the  actual  loading  distribution.  However,  since  the  force  records  were 
used  to  determine  that  load  distribution,  it  is  subject  to  experimental 
error. 

The  discrepancy  between  the  theoretical  and  experimental  results,  as 
well  as  the  scatter  in  the  experimental  results,  may  be  due  to  the  following 
considerations : 

1)  The  nondimensional  buckling  pressure  is  a  function  of  the  aspect 
ratio  (structure  width/characteristic  length).  Thus,  the  scatter  in  the  ex¬ 
perimental  data  can  be  partially  attributed  to  the  fact  that  the  aspect  rati¬ 
os  were  not  the  same  for  each  test.  The  discrepancy  between  the  theoretical 
and  experimental  results  can  be  partially  explained  by  the  fact  that  the 
aspect  ratio  used  in  the  theoretical  analysis  was  rarely  equal  to  the  aspect 
ratios  measured  during  the  experimental  tests.  Since  a  10%  error  in  the 
characteristic  length  (the  structure  width  was  a  constant)  results  in  a  20% 
error  in  the  nondimensional  buckling  pressure,  the  differing  values  used  for 
the  characteristic  lengths  can  be  a  significant  source  of  error. 

2)  In  the  theoretical  analysis  the  ice  sheet  is  assumed  to  be  a  perfect 
structure;  i.e.  there  are  no  imperfections  such  as  initial  curvatures,  eccen¬ 
tricity  of  loads,  etc.,  in  the  ice  sheet.  Thus,  the  buckling  loads  computed 
in  the  theoretical  analysis  are  the  loads  at  which  bifurcation  of  equilibrium 
exists  for  "perfect"  ice  sheets.  Any  imperfection  in  the  sheet  or  eccentri¬ 
city  in  the  loading  can  lead  to  lower  buckling  loads  in  the  experiments. 

The  authors  believe  that  the  scatter  in  the  experimental  data  can  large¬ 
ly  be  attributed  to  the  variations  in  the  ratio  of  structure  width  to  charac- 
teristic  length,  as  the  buckling  load  (P/BKL  )  is  most  sensitive  to  this  par¬ 
ameter  (Fig.  15). 
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